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Analysis of homogeneous turbulent reacting flows 

By A. D. Leonard 1 , J. C. Hill 1 , S. Mahalingam 2 , AND J. H. Ferziger 2 


Full turbulence simulations at low Reynolds numbers were made for the single- 
step, irreversible, bimolecular reaction between non-premixed reactants in isochoric, 
decaying, homogeneous turbulence. Various initial conditions for the scalar field 
were used in the simulations to control the initial scalar dissipation length scale, 
and simulations were also made for temperature-dependent reaction rates and for 
non-stoichiometric and unequal diffusivity conditions. Joint probability density 
functions (pdf’s), conditional pdf’s, and various statistical quantities appearing in 
the moment equations were computed. Preliminary analysis of the results indicates 
that compressive strain-rate correlates better than other dynamical quantities with 
local reaction rate, and the locations of peak reaction rates seem to be insensitive 
to the scalar field initial conditions. 


1. Background 

Chemically reacting species that are not premixed are separated by a reaction 
zone, unless the time scale for mixing is much less than the time scale for reaction. 
The structure of the reaction zone is determined by interactions of the turbulent 
motion, molecular diffusion, and chemical kinetics. Molecular diffusion controls the 
rate of reaction when the kinetic rate is very fast and the reaction zone becomes a 
surface separating regions of reactants. This surface is typically a contour surface 
for a conserved scalar species. The reaction zone has a finite thickness for fast but 
finite rate kinetics. The reactants become completely mixed when the kinetic rate 
is very slow, and the classification of “non-premixed” is then irrelevant, because the 
kinetics, rather than transport processes, control the rate of reaction. 

Previous simulations of inert scalars are useful in the study of reacting flows 
because an infinitely fast reaction can be described in terms of a conserved scalar 
species. Simulations of inert scalars in homogeneous turbulence show a tendency for 
the scalar gradient to align with the most compressive rate of strain of the velocity 
field and for the rate of scalar dissipation to be highest where the alignment of the 
scalar gradient and the direction of the compressive strain rate is greatest (Ashurst 
et aL, 1987). Previous work with reacting scalars (Leonard and Hill 1987, 1988; 
Hill, Leonard & Rogers 1987) show that the compressive strain rate direction tends 
to be perpendicular to the reaction zone, in agreement with the results of Ashurst 
et at. 
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Probability density {unctions (pdf’s) are frequently used in the statistical treat- 
ment of reacting turbulent flows, because the source terms do not need to be directly 
modeled. Mixing effects are seen in the statistical equations as conditional expec- 
tations of scalar dissipation. Eswaran and Pope (1988) have used direct numerical 
simulation to study the evolution of the pdf and the conditional scalar dissipation 
of an inert scalar species. The initial mixing rate was increased by a decrease in 
length scale, but evolution of the shape of the pdf of the scalar was not affected. 
The pdf was initially a “double delta” function, which would be appropriate for the 
conserved scalar of unmixed reactants. It evolved into a nearly Gaussian form. The 
conditional dissipation approached a self-similar value and became independent of 
the level of the scalar. The conditional dissipation needs to be modeled properly to 
give accurate mixing results, but current modeling efforts ( e.g Pope, 1985; Koll- 
mann, 1987) do not account for the effects of reaction rate on the mixing terms. 
The price to pay for treating the source terms exactly in reacting flows is the need 
to model all of the mixing effects, and this must be done in a higher dimensional 
system (Hill, 1988). 

2. Problem Description 

The specific problem under study is an irreversible reaction of two initially segre- 
gated species in homogeneous turbulence. The reaction mechanism can be expressed 
as 

A + B — ► P. 

Full turbulence simulations were performed for the evolution of the concentration 
and velocity fields, taking into account the advection and Fickian diffusion of species 
A, B, and P. The turbulence in most cases was isotropic and decaying. The velocity 
field was incompressible and not affected by the concentration fields, even in the 
studies using a temperature-dependent reaction rate. The governing equations for 
the study are 

8A 

— + u* VA = P a Vo2A + uu, 

ot 

with similar equations for the evolution of species B and P, and 

d\i Vp _ 

— + u • Vu = + i/V o 2u. 

ot p 

The concentration of species A is denoted by A and the rate of change due to 
reaction by wa • Two forms of the source terms were used in the studies, 

wa — — kAB 

w A = -k 0 eo-T/T a AB, 

where the k' s are reaction rate constants and T a is the activation energy in the 
nonisothermal case. 

The initial concentration distribution was one of two forms: alternating stripes 
of reactant species aligned perpendicular to one of the coordinate axes, or spatially 
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segregated species uniform on the largest scales of the domain. The second form 
of initial conditions was created using a procedure similar to that used by Eswaran 
and Pope (1988). The initial velocity field was chosen by randomly selecting Fourier 
coefficients and scaling them so the three dimensional energy spectrum fit a specified 
function. This initial velocity field was allowed to develop in a “pre- simulation” until 
the skewness of the velocity derivative had a value that agreed with experimentally 
measured conditions in turbulent flows. 

The simulation method used an existing pseudospectral code (Rogallo, 1981) that 
was modified to include the rate of change of scalar values due to chemical reactions. 
The domain for most studies was 64o3 points. Three velocity components and eight 
scalar values were included in the numerical solutions. 

One objective of the study was to determine the influence of local properties 
of the turbulence on the local reaction rate and structure of the reaction zone. 
This includes studies of the alignment of directions of strain rate and concentration 
gradients for reacting scalars for different speeds of reaction and of the influence of 
vorticity on reaction zone structure. The long-term goal of such a study is to model 
average reaction rate using strain-rate parameters. 

A second objective was to determine the importance of the initial conditions in 
previous calculations (Leonard and Hill 1987). In these studies with a preferential 
alignment of the scalar fields and no variation in scalar length scales, globally- 
averaged scalar microscales showed little sensitivity to reaction rate. 

A third objective was to determine if there is any effect of reaction rate on con- 
ditional expectations of molecular diffusion terms in a hierarchy of pdf equations. 
The modeling of these terms neglects any dependence on reaction rate, despite the 
fact that the pdf approach is widely used in the statistical treatment of reacting 
flows. 

Additional questions concerned the behavior of reactions with non- stoichiometric 
conditions or with species having different mass diffusivities, as well as the local 
structure of the reaction zone for temperature-dependent kinetics. 

2. Approach 

The existing databases of homogeneous turbulent flows with chemical reaction 
were inadequate to answer the questions that had been raised, and so the emphasis 
was placed on the generation of new databases during the Summer Program, with 
careful consideration of physical parameters and initial conditions. Additional diag- 
nostics were developed and added to the code to permit evaluation of the evolution 
of single point pdf’s, single-point conditional expectations, and averages over two, 
rather than three, spatial dimensions during the course of the simulations. 

The conditions of the studies made during the Summer Program are summarized 
in Tables I and II. Table I groups the simulations into 7 different studies, each 
study consisting of one or more simulations. Parameters for the scalar variables 
were varied within each study, as well as between different studies, in order to 
isolate effects such as initial conditions or kinetic mechanism. 

Most studies have in common a decaying isotropic turbulent velocity field which 
has evolved from what will be called a “developed” initial field, isothermal reaction, 
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Table I. Summary of Studiesoa 


Name of Study Effects Studied Initial Conditions Comments 

1. “Non-ideal” 

a) Non-stoichiometric 

b) Unequal diffusivity 

Migration of 
reaction zone 

Gaussian v-field 
(Same as ’87 study) 
Scalar slabs 

Reaction zone 
follows X = 0 

2. “Slab” 

Kinematics of reaction 
zone (Improved v-field) 

Developed v-field 
Scalar slabs 

Confirmed ’87 results 
re: strain effects 

3. “Flame” 

Ignition and 
quenching 

Developed v-field 
Scalar slabs 
Isothermal 

Arrhenius rate expression 
Quenching not observed 

4. “Stirred” 

Effect of scalar field 
initial conditions 

Developed v-field 
Prestirred scalar 

Conditional scalar 
dissipation measured 

5. “Stripes” 

Effect of scalar 
dissipation scale 

Developed v-field 
Scalar stripes 
(1,2,3, 4) 

Conditional scalar 
dissipation measured 

e. “Forced” 

Decouple reaction rate 
from turbulence decay 

Developed v-field 
Scalar slabs or other 

Incomplete (Compiled 
code untested) 

7. “Shear” 

Strain-rate dependence 
on reaction rate and 
scalar field structure 

Top-hat v-spectrum 
Scalar slab in 
plane of shear 

Incomplete (Excessive 
Gibbs’ ringing for 
conditions used) 


Oa All studies were for decaying homogeneous turbulence except for #6 and 7. The term “scalar 
slabs* denotes one stripe of one reactant centered in the domain of non- premixed reactants, 
hence there are two reaction zones initially. Improved diagnostics, including various joint-pdf’s 
and planar averages, were used for all studies except #1. 


equal mass diffusivities for all species, and stoichiometric proportions of reactant 
concentrations. The parameters and initial conditions shared by a majority of the 
simulations are given in Table II. The individual features of each study that were 
changed are discussed in the remainder of this section, using the informal names in 
Table I as convenient identifiers. 

1. “Non- ideal.” The non-ideal study examined the migration of the reaction zone 
for two different conditions. One simulation was made with stoichiometric coeffi- 
cients (Bq/Aq) of 1 and 2 for two different reaction rates. The second simulation 
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Physical parameters 

V 

0.02 

Vo a 

0.0286 

k } ko 

variousO b 

Initial valuesO c 

u' 

1.03 

^ 9 

0.381 

A/ 

0.844 

A,B 

1.0 

Xa 

various 


Oa Da = Dq = T) for all studies except #lb. 

ob The parameter was varied within some studies, as well as be- 
tween studies. 

Oc These initial conditions are for the “developed” velocity field 
used in studies #2-5, and for a stoichiometric proportion of 
reactant concentrations. 

used a mass diffusivity ratio, ( Da/Db )> of 1 and 3.5 for a stoichiometric ratio of 
1. The initial reactant concentrations were the “scalar slabs,” meaning one stripe 
of each reactant. The initial energy spectrum for the velocity field fit a Gaussian 
form, and the velocity field was allowed to develop for 100 time steps before the 
reaction began. The initial value of R\ was about 20. The Damkohler numbers 
based on the diffusivity and initial mean concentration of species A and defined 

by 

Dai = kAoAj/u 1 

Dan = IcAqXa o 2/Da 

were 0.8 and 4 for the first kind and 30 and 150 for the second kind. 

“Slab.” The slab study was aimed at repeating the results from the previous 
Summer Program with a different detailed velocity field. Fourier coefficients for 
the velocity field were chosen from a Gaussian distribution, subject to satisfying 
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the continuity equation and an initial energy spectrum E(k) oc k o4/(l + ( kL ) o 
3)o2, and allowed to develop for 100 time steps. This initial field had more energy 
at high wavenumbers than did the field used for the 1987 study, and so R\ was 
smaller. The initial value of R\ was 19.7. Two reactions were included in the 
simulation with different reaction rate coefficients. The Damkohler numbers for 
this study were 1.64 and 6.56 for the first kind and 62.8 and 251 for the second 
kind. Eight scalars were carried in the simulation: the reactant and product 
concentrations for each reaction and the concentrations of two inert species having 
different initial conditions. 

3. “Flame.” The “flame” study used the same initial conditions for the velocity and 
reactant species as in the slab case, but the reaction rate was temperature depen- 
dent, with an Arrhenius-type rate constant. For this case the temperature was 
calculated (actually, the fractional approach to the adiabatic flame temperature 
was used), rather than the product concentration. Simulations were performed 
for two different values of the Lewis number, in order to examine the possibility 
of strain-induced quenching, and also for two different values of the Zeldovich 
number. The Zeldovich number is defined here as the logarithm of the ratio of 
the reaction rate at the adiabatic flame temperature to the reaction rate at the 
initial temperature. Zeldovich numbers of 1 and 2 were used in the study. The 
Damkohler numbers based on the initial temperature and mean concentrations 
were the same as those used in the “slab” study. 

4. “Stirred*” Conditions for the “stirred” study were identical to the “slab” case, 
except for having different initial scalar distributions. A scalar field was defined in 
a manner similar to that for the velocity field, by scaling Fourier coefficients to fit 
an energy spectrum and advancing the coefficients in a “pre- simulation”. In the 
case of the reactant concentrations, one reactant was set to be a (positive) initial 
value at the points were the initial scalar was positive and the other reactant 
was set to be a (positive) initial value at the points where the initial scalar was 
negative. The high wavenumbers were damped to eliminate Gibbs’ ringing in the 
initial conditions. The reactant species were spatially segregated, but uniform 
on a large scale. The microscales for the reactants were much smaller than in 
the slab case, so mixing will tend to be faster. The Damkohler numbers of the 
second kind were 7.71 and 30.9, while the Damkohler numbers of the first kind 
were identical. 

5. “Stripes.” The initial values of reactant concentration in the slab and stirred 
cases differed both in orientation of the reaction zones and in the dissipation 
length-scale. The stripe study isolates the effect of the initial concentration scale 
from the orientation. One reaction rate coefficient is used for the calculation of 
reactant species for 4 reactions, each with different initial conditions. Damkohler 
numbers of the second kind were 13.7, 19.2, 30.4, and 62.8. 

6. “Forced*” Reactions in decaying homogeneous turbulence are unduly influenced 
by the amount of mixing in the early part of the simulation and controlled by 
diffusion in the later parts of the simulation. A forcing mechanism was devised 
to maintain a constant energy level but has not been implemented at this point. 
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7. “Shear.” The same scalar field initial conditions of the “slab” case were used 
in the shear case, with the initial reaction zones lying in the plane of shear, 
but with a mean velocity gradient imposed on the flow. The initial velocity 
field was defined using the top-hat energy spectrum E(k) = constant, for k t < 
k < & 2 - The resolution was doubled in the direction of the mean flow, using 
a wavenumber grid 64 x 64 x 128. A preliminary run was made — at too high 
a Damkohler number, apparently — but discontinued because of excessive Gibbs’ 
ringing. Intensification of vorticity that occurs in the shear case clearly imposes 
limitations on the Damkohler number range that can be simulated. 

3. Preliminary Results 

Because most of the effort at the Workshop was spent in developing new databases 
for subsequent analysis, very little analysis was accomplished at the Workshop itself, 
and the preliminary results discussed here were obtained by cursory examination 
of some representative save-fields. Interpretation of the statistical results, espe- 
cially the dependence of conditional expectations of scalar dissipation rate on the 
Damkohler number (dimensionless reaction rate constant), is still in progress, as is 
a more detailed examination of the saved fields. The principal results obtained so 
far deal with the kinematics of the reaction zone and with the influence of initial 
concentration length scales on the reaction (Hill & Leonard, 1988; Leonard et ai, 
1988), and some of those results are summarized below. 

The main observation from the “Non-ideal” study was the movement of the 
reaction zone for both the non-stoichiometric and the unequal diffusivity cases. 
The maximum reaction rate for the non-stoichiometric case occurred at the same 
point as the value of the conserved scalar corresponding to stoichiometric conditions 
(X = A — B = 0). The conserved scalar treatment is not valid for the case of 
reactants with unequal diffusivities, but the behavior is qualitatively the same for 
the two cases. The mass flux of reactants to the reaction zone must be the same 
for both reactants, but the fluxes are unequal in the non-ideal cases, either because 
of a difference in diffusivity or in concentration gradient, and so the reaction front 
must move. 

“Slab” initial conditions are used to study the structure of the reaction zone 
because two distinct fronts are present in the simulation. In the present study 
the initial reaction zones are perpendicular to the y axis and located at y = 7r/2 
and y = 37r/2. The average positions of the fronts will not change in isotropic 
turbulence for reactions between stoichiometric proportions of reactants with equal 
diffusivities. 

The scalar field has a preferential alignment with this initial condition. The gradi- 
ents of reactant concentrations are in the direction of the y axis. The initial velocity 
condition is isotropic, so the probability of the cosine of the angle between concen- 
tration gradient and any of the eigenvectors of the rate of strain tensor is uniformly 
distributed. The concentration gradients become aligned with the direction of most 
compressive strain. Figure 1 shows the pdf’s of the angles between the gradient of 
concentration for species A and eigenvectors of the rate of strain tensor for the slab 
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FIGURE 1. Probability density function /(cos0) of the alignment of the gradient 
of reactant species concentration and the eigenvalues of the strain-rate tensor, for 
a slab case at t — 2.75. 


case at a dimensionless time of 2.75. The pdf is sharply peaked at cos 0 7 = 1 (where 

is the angle between the gradient of species A and the eigenvector corresponding 
to the most compressive strain rate) after only 100 steps or t — 0.4. The alignment 
is in agreement with the results of Ashurst et al. (1987), but the present results 
show how rapidly it develops. 

A cross section of the domain in the slab case, either in an x—y or a y—z plane, 
shows two distinct, distorted reaction zones, as seen in Figure 2 . The degree of 
alignment with the eigenvector corresponding to 7 , the least principal strain rate, is 
shown in Figure 3. Concentration gradient alignment with compressive strain is im- 
portant because the strain increases the magnitude of the gradient, which increases 
the mass transport of reactant species to the reaction zone, which increases rate of 
reaction. This is shown by the set of equations developed from the conservation 
equations for the reactants in the case of an isothermal reaction, 


dA 

dt 

dAB 

dt 


= -k AB, 


= -2VVA . VB - k AB(A + B ), 


dVA ' — = -2 VVVA : VVfl - 2VA • e • VJ 3 

dt 

- k [(A + B)VA • VB + BVA • VA + AVB • VB . 
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We refer to the term VA*e* VB as “gradient compression”, because the gradients are 
being amplified by the tendency of turbulence to move material surfaces together. 
When the scalar gradients are aligned with 7 the gradient compression becomes 
— 7 |VA||VJ?|. This term, which is positive because the gradients of A and B are 
in opposite directions and 7 < 0, increases the magnitude of V A • VJ5 at a rate 
e o —2yt in the absence of diffusion or reaction. Contours of AB, |VA * VjB|, and 
VA-e-VB all show peak values in the same regions. 

A second consequence of the scalar gradient alignment is that the location of peak 
reaction rates seems not to be very sensitive to the initial conditions of the scalar 
field. The initial values of the ‘‘Stirred” case were determined randomly, while the 
initial values of the “slab” case were chosen deterministically. The velocity field was 
the same for the two cases. A cross section of the slab case in a x—z plane, as well 
as cross sections of the stirred case in any plane, show a more random structure 
than that in Figure 2. Regions of intense reaction rate for the slab and stirred cases 
coincide in a plane that was the center of one of the initial reaction zones for the 
slab case (Figs. 4 and 5), even though the initial concentration fields are different. 

Since the reaction rate is high in the same locations for completely different initial 
conditions, it may be possible to predict the structure of the reaction zones from a 
knowledge of velocity field. One parameter that must be used is the rate of strain. 
Criteria were proposed by J. Hunt and A. Wray at the CTR Summer Program to 
determine kinematic structures of the flow field. These structures are identified as 
streams, eddies, and convergence zones, based on the strain rate invariants, vorticity, 
kinetic energy, and pressure. Hunt and Wray suggest that the reaction rate will be 
highest in the convergence zones. Simulation results, however, show that regions 
where the reaction rate is high are not necessarily convergence zones. The gradient 
compression term appears to be the best marker of regions of high reaction rate. 

Only a few remarks will be made here about the remaining studies. First, the 
“Flame” study. The reaction rate has been found to be higher when the magnitude 
of 7 is large, because of the gradient compression. Gradient compression is believed 
to be responsible for quenching in flames, because the transport of heat away from 
the reaction zone becomes more important than the transport of reactants to the 
reaction zone for Arrhenius kinetics. Simulations with temperature-dependent re- 
action rate do not show any reduction of reaction rate where the strain rates are 
high. The only qualitative difference between the temperature-dependent case and 
the isothermal case is that the reaction zones are initially thinner in the former 
case, an apparent ignition effect. 

Comparison of results from the “Stripes” study with those of the “stirred” 
and “slab” studies shows that as one increases the scalar dissipation scale, mixing 
becomes more important, and the reactants decay faster. (The fraction of the 
domain occupied by scalar dissipation zones, and hence reaction zones, is increased.) 
In the case of several stripes, interference effects clearly become important after 
several hundred time steps, and the reaction zones lose their identity. 

The studies of chemical reactions in “Forced” and “Shear” turbulence were 
not carried out during the Workshop (except for coding of the forcing algorithm in 
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FIGURE 4. Contours of reaction rate, JcAB, for a slab case at t = 2.75 in the plane 
of the original reaction zone, y = 7r/2. These data are for the same time and the 
same run as figure 3; only the planes observed are different. 



FIGURE 5. Contours of reaction rate, kAB } for a stirred case at t = 3.12 in the 
plane y = 7r/2. Coincidence, at nearly the same time, of the locations of the highest 
reaction rates as for the slab case in figure 4 should be noticed. 
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the former case and making a preliminary run at too high a Damkohler number in 
the latter), and those studies will be continued. 

4. Conclusions 

The local rate of strain in the fluid appears to be the feature of turbulence that has 
the most direct influence on the course of a chemical reaction, indicating that rate 
of strain should be incorporated into models for reacting flows. Reaction zones are 
aligned perpendicular to the direction of compressive strain rate, and the magnitude 
of the reaction rate is greatest when the magnitude of the strain rate is highest. 
The straining motion of the fluid compresses concentration gradients to enhance 
molecular diffusion, and thereby reaction rate. These observations support the 
simple one-dimensional model by Gibson and Libby (1972) of the reaction zone, in 
which the reaction surface is aligned with the most compressive strain rate, and 
the fluid motion is approximated by u= — jx in a coordinate system fixed to the 
reaction zone. 

Preliminary examination of the saved fields suggests furthermore that the align- 
ment of scalar gradient and strain rate directions is independent of initial conditions, 
and that the highest reaction rates tend to occur at the same points in the flow for 
reactions with different initial reactant distributions. Changing the initial scalar 
microscale, and hence changing the fraction of the domain occupied by reaction 
zones, affects the overall reaction rate and scalar dissipation rates. 

The analysis of the data generated during the Summer Program is continuing, 
and additional simulations are being made where necessary. Further studies of 
non-isothermal cases are planned, as well as non-decaying turbulent flows. 
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